function J = analJacobian(Z,x_t,y_tN,N,setupEPer)
% Recall we need the Jacobiant for vecZ where
% vecZ = [y_t;
%
%         x1_t+1;
%         y_t+1;        
%
%         x1_t+2;
%         y_t+2;;
%         
%         ...
%
%         x1_t+N-1
%         y_t+N-1
%
%         x1_t+N]

% Unfold setupEPer
nx       = setupEPer.nx;
mx       = setupEPer.mx;
myx      = setupEPer.myx;
ny       = setupEPer.ny;
params   = setupEPer.params;
h0       = setupEPer.h0;
hx       = setupEPer.hx;
g0       = setupEPer.g0;

[x,y]    = VarStack(Z,nx,mx,myx,ny,hx,x_t,y_tN,N);
J        = zeros((ny+mx)*N,(ny+mx)*N);
getDeriv = 1;

for t=1:N 
    [~,fx1,fx1p,fy,fyp,fx12] = getModelDeriv_EP(nx,mx,myx,ny,params,N,h0,g0,x,y,t,getDeriv);
    
    if t == 1
        J(1+(ny+mx)*(t-1):(ny+mx)*t,1:2*ny+mx) = [fy fx1p fyp];
    elseif t < N
        J(1+(ny+mx)*(t-1):(ny+mx)*t,(ny+mx)*(t-2)+1:(ny+mx)*(t-2)+myx) = fx12old;
        J(1+(ny+mx)*(t-1):(ny+mx)*t,ny+(ny+mx)*(t-2)+1:ny+(ny+mx)*t)   = [fx1 fy fx1p fyp];
    elseif t == N
        J(1+(ny+mx)*(t-1):(ny+mx)*t,(ny+mx)*(t-2)+1:(ny+mx)*(t-2)+myx)      = fx12old;
        J(1+(ny+mx)*(t-1):(ny+mx)*t,ny+(ny+mx)*(t-2)+1:ny+(ny+mx)*(t-1)+mx) = [fx1 fy fx1p]; 
    end
    fx12old = fx12;
end
J = sparse(J);
%length(nonzeros(J))
end


